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ABSTRACT 

This paper continues our previous exploration of the effects of turbulence on mean 
motion resonances in extrasolar planetary systems. Turbulence is expected to be present 
in the circumstellar disks that give rise to planets, and these fluctuations act to com- 
promise resonant configurations. This paper extends previous work by considering how 
interactions between the planets and possible damping effects imposed by the disk af- 
fect the outcomes. These physical processes are studied using three related approaches: 
direct numerical integrations of the 3-body problem with additional forcing due to tur- 
bulence, model equations that reduce the problem to stochastically driven oscillators, 
and Fokker-Planck equations that describe the time evolution of an ensemble of such 
systems. With this combined approach, we elucidate the basic physics of how turbulence 
can remove extrasolar planetary systems from mean motion resonance. As expected, 
systems with sufficiently large damping (dissipation) can maintain resonance, in spite of 
turbulent forcing. In the absence of strong damping, ensembles of these systems exhibit 
two regimes of behavior, where the fraction of the bound states decreases as a power-law 
or as an exponential. Both types of behavior can be understood through the model de- 
veloped herein. For systems that have weak interactions between the planets, the model 
reduces to that of a stochastic pendulum, and the fraction of bound states decreases as a 
power-law Pb oc For highly interactive systems, however, the dynamics are more 

complicated and the fraction of bound states decreases exponentially with time. We 
show how planetary interactions lead to drift terms in the Fokker-Planck equation and 
account for this exponential behavior. In addition to clarifying the physical processes 
involved, this paper strengthens our original finding that turbulence implies that mean 
motions resonances should be rare. 



Subject headings: MHD — planetary systems — planetary systems: formation — plan- 
ets and satellites: formation — turbulence 



^Michigan Center for Theoretical Physics 
Physics Department, University of Michigan, Ann Arbor, MI 48109 

^Physics Department, University of Wisconsin, Madison, WI 53706 

^Astronomy Department, University of Michigan, Ann Arbor, MI 48109 

^Department of Mathematics, University of Michigan, Ann Arbor, MI 48109 



- 2 - 



1. INTRODUCTION 

A sizable fraction of the observed extrasolar planetary systems with multiple planets have 
period ratios that are close to the ratios of small integers, and hence are candidates for being in 
mean motion resonance (e.g., Mayor et al. 2001, Marcy et al. 2002, Butler et al. 2006). However, 
a true resonance not only requires particular period ratios, but also oscillatory behavior of the 
resonant angles, which in turn requires specific values for the orbital elements of the two planets 
(Murray & Dermott 1999; hereafter MD99). As a result, a mean motion resonance represents a 
rather special dynamical state of the system. The suspected origin of such configurations is through 
a process of convergent migration, wherein two planets are moved inward together through their 
solar system by a circumstellar disk and thereby given the opportunity to enter into a resonant 
state (e.g., Lee & Pealc 2002, Lee 2004, Bcaugc ct al. 2006, Moorhead & Adams 2005, Crida et al. 
2008). Since these disks are likely to be turbulent (e.g., Balbus &: Hawley 1991), and since resonant 
states are easily disrupted, it is possible for turbulent fluctuations to drive planetary systems out 
of mean motion resonance. In an earlier paper (Adams et al. 2008; hereafter Paper I), we explored 
this possibility and found that the presence of turbulence implies that mean motion resonances 
should be rare. This paper continues this earlier effort by exploring the problem in greater depth 
and by including additional physical effects. 

This work is motivated by the current observations of extrasolar planetary systems with multi- 
ple planets. Although many of the observed systems display period ratios that are consistent with 
mean motion resonance, more data is required to determine if many of these systems are truly in 
resonance. The review of Udry et al. (2007) shows that four systems have 2:1 period ratios, 2 
systems have 3:1 period ratios, 2 systems have 4:1 period ratios, and 4 systems have 5:1 period 
ratios. Of these possible resonant systems, those with 2:1 period ratios are the most compelling 
candidates, with ratios P2/P1 = 2 it 0.01. Of these systems, one case is GJ876, which is well 
observed and has two planets thought to be deep in a 2:1 mean motion resonance (Marcy et al. 
2001, Laughlin &: Chambers 2001). The other candidates for systems in 2:1 resonance include HD 
82943 (Gozdziewski & Maciejewski 2001, Mayor et al. 2004, Lee et al. 2006), HD 128311 (Vogt 
et al. 2005), and HD 73526 (Tinney et al. 2006, Sandor et al. 2007). The status of these latter 
three systems remain unresolved; in any case, their libration widths are wide (as well as uncertain) 
and hence these systems are not as "deep" in resonance as the GJ876 system. In addition, the 
55 Cancri system (Marcy et al. 2002) has planets with period ratios close to 3:1, and has been 
suggested as another candidate for mean motion resonance (Ji ct al. 2003); however, subsequent 
work (Fischer et al. 2008) indicates that the system is unlikely to be in resonance. Finally, we 
note that a new system of "Super-Earths" (with masses of 4.2, 6.9 and 9.2 Me) has recently been 
discovered orbiting HD 40307 with periods of 4.3, 9.6, and 20.5 days (Mayor et al. 2008). This 
system thus has period ratios that are relatively close to 4:2:1, but are most likely too distant 
to be in resonance. Although more data is required, the observational landscape can be roughly 
summarized as follows: one system (GJ876) is deep in a 2:1 mean motion resonance, three more 
systems are either in resonance or relatively close, and many more planetary systems have period 
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ratios indicative of mean motion resonance but are probably not actually bound in resonant states. 
Since mean motion resonances are relatively easy to compromise, this collection of observed systems 
provides important clues regarding their formation and past evolution. 

In Paper I, we considered turbulent fluctuations as a mechanism to remove systems from mean 
motion resonances while retaining nearly-integer period ratios. If the fluctuations have sufficiently 
large amplitude and duty cycle, turbulence is indeed effective at removing systems from bound 
resonant states. Specifically, this earlier work used both a simple model equation (a stochastic 
pendulum) and direct integrations based on the architecture of the observed planetary system 
GJ876 (see above). In both cases, turbulent forcing, with the amplitudes expected from magneto- 
hydrodynamic (MHD) simulations (e.g.. Nelson & Papaloizou 2003; Laughlin et al. 2004, hereafter 
LSA04; Nelson 2005) of magneto-rotational instability (MRI), were found to drive systems out of 
a resonant state. We note that not all circumstellar disks will support MRI, especially when they 
are sufHciently resistive. Nonetheless, this treatment applies to any type of stochastic fluctuations, 
including all types of turbulence, that could be present in such disks. Although our earlier work 
elucidates the basic physics of the problem, and shows that turbulence can be effective at removing 
systems from resonance (Paper I), several issues must be studied further: 

[1] The planets migrate inward during much of the time while turbulence is expected to be 
present. As a result, the effects of turbulence on mean motion resonance during the migration 
epoch must be considered. This issue has been considered in recent work (Moorhead 2008), which 
shows that turbulence that acts during the epoch of planet migration does indeed compromise mean 
motion resonance in a manner that is qualitatively similar to the results of Paper I; a more detailed 
study of this process is underway (Moorhead &; Adams, in preparation). 

[2] The circumstellar disk that produces both the inward migration and the turbulent fluctu- 
ations can also produce a net damping of the eccentricity of the planetary orbits. This damping, 
in turn, can act as a damping mechanism to counter the action of turbulence driving systems from 
resonance. The process of convergent migration itself, which allows planets to enter into mean mo- 
tion resonance (e.g., Lee & Peale 2002; Crida et al. 2008), can also provide a damping mechanism 
for the resonance. More generally, damping can effectively take place whenever the fluctuations 
have a nonzero flrst moment (R. Malhotra, private communication) or through dynamical viscosity 
(E. Zweibel, private communication). In any case, the combined effects of damping and turbulent 
fluctuations on mean motion resonance should be studied further, and are considered here (see §3). 

[3] The stochastic pendulum model (used in Paper I) reduces the full physical system — which 
has 18 phase space variables for a two planet solar system — to a much simpler system with only 
one variable. In its simplified form, the system can freely random walk both in and out of resonance 
through the action of the turbulent fluctuations. As mentioned above, however, a true resonance 
represents a special dynamical state of the system, and hence the full system (with 18 phase space 
variables) will not necessarily re-enter resonance as easily as the one variable system. In Paper 
I, we found that the full system can in fact random walk back into resonance after leaving, but 
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not as easily as suggested by the stochastic pendulum model. As a result, the issue of how a 
planetary system re-enters into resonance after leaving, especially in the presence of fluctuations, 
remains open. Any barrier to rc-cntry into resonance causes the fraction of bound states (systems 
in resonance) to decrease with time faster than the stochastic pendulum model (see §4). A related 
issue is that highly interactive systems, those with large planets and high eccentricities, tend to 
leave resonance more readily than systems that are less interactive. 

[4] When planets are in resonance, they are protected (at least in part) from orbit crossings, 
even when the eccentricities are fairly large. When the planets are knocked out of resonance, 
however, they can be subject to greater eccentricity excitation and eventual orbit crossing, and 
hence have a chance to scatter off each other and be ejected. As a result, solar systems that leave 
resonance have a greater chance of losing planets. Further, when a system ejects a planet, that 
system can never re-enter a resonant state; this loss of planets contributes to the decreasing number 
of bound states (systems in resonance) as a function of time. We also expect highly interactive 
systems to eject planets more readily than those that are less interactive. The role of planet ejection 
on the time evolution of the fraction of bound states must be explored further. 

The goal of this paper is to study the effects of turbulence on mean motion resonance with 
the above issues in mind. In §2, we review the three approaches to the problem used herein: The 
basic model of the resonant system as a (stochastic) pendulum, the corresponding phase space 
approach to the dynamics using a Fokkcr-Planck equation, and full numerical integrations (using 
all 18 phase space variables for the 3-body problem). This section also presents new ensembles of 
numerical simulations to illustrate how well the systems are described by the stochastic pendulum 
model and where its shortcomings are found. In §3 we consider the combined effects of damping 
and turbulence on mean motion resonance. In §4 we derive a more detailed model of mean motion 
resonance, where we keep higher order terms that include effects from planet-planet interactions. 
Finally, we conclude in §5 with a summary of our results and a discussion of their implications. 

2. FORMULATION 

In this section we outline the three basic methods used in this paper. We first review the sim- 
plest version of the pendulum model for mean motion resonance (from MD99; see §2.1) and then 
outline our treatment for including turbulent fluctuations (from Paper I; see §2.2). Our first ap- 
proach is to directly integrate the resulting stochastic differential equations that describe the time 
evolution of resonances in the presence of turbulence; in this case, a large ensemble of different 
solutions must be explored to sample the effects of different realizations of the stochastic fluctua- 
tions. An alternate approach, outlined in §2.3, is to find the corresponding Fokker-Planck equation, 
which describes the time evolution of the distribution of states. Our third and final technique is to 
directly integrate (numerically) the full 3-body system including velocity perturbations to account 
for the turbulence (§2.4); this approach also requires a large ensemble of different realizations of 
the problem to build up a statistical description. 
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2.1. The Basic Pendulum Model 

As a starting point, we use a version of the formalism presented in MD99. To start, this 
treatment is limited to the circular restricted 3-body problem and hence considers the libration 
angle (f) for two planets in resonance to have the form 

= jlA2 + j2A + j4ro, (1) 

where A and A2 arc the mean longitudes of the two planets and w is the longitude of periapse of 
the inner planet. In this case, the equation of motion for the resonance angle (f) reduces to that of 
a pendulum, i.e., 

^+u;2sm0 = O. (2) 

The natural oscillation frequency loq of the pendulum — the libration frequency of the resonance 
angle — is given by 

= -Sj^Crne^^^^ , where Cr = nVtafd{oi) . (3) 

Here, O and e are the mean motion and eccentricity of the inner planet, and (j2, J4) are integers 
that depend on the type of resonance being considered. The parameters Cr, e, and f2 are assumed 
to be constant for purposes of determining the frequency. The mass ratio ji = m2/M^, where m2 
is the mass of the outer planet and is the stellar mass. The quantity afd[a) results from the 
expansion of the disturbing function (MD99), and the parameter a is the ratio of the scmimajor 
axes of the two planets, i.e., a = ai/0,2. Note that this approximation scheme (adapted from MD99) 
neglects terms of order 0{fi); we consider higher order terms in §4. 

Many of the observed (candidate) resonant systems are in or near the 2:1 mean motion res- 
onance; for the sake of definiteness, we focus on that case. This resonance is also generally the 
strongest. For the 2:1 mean motion resonance, the integer j2 = —1 and afd{a.) —3/4 (MD99). 
In this case, the natural oscillation frequency loq of the libration angle is given by 

c.2«^^el^4lo2. (4) 

Typical planet masses in observed extrasolar planetary systems are of order a Jovian mass so that 
~ 10~^. The eccentricity can vary over a wide range, with e ~ 0.1 in order of magnitude, but the 
median eccentricity is closer to e ~ 0.3. The relevant values of |j4| = (1,2), and we generally use 
= 2 (see the discussion of §2.4). As a result, the ratio of frequencies ujq/Q, ~ 10^^ and hence 
the period of the libration angle is expected to be ~ 100 orbits. With this frequency specified, we 
define a dimensionless time variable 

T = LOot, (5) 

and write the pendulum equation in dimensionless form 

^+sin(^ = 0, (6) 
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Note that the energy scale associated with the potential well of the resonance, modeled here 
as a pendulum, is much smaller than the binding energy of the planets in their orbits. Neglecting 
dimensionless numbers of order unity, we can write the specific energy Ep of the pendulum as 

Ep ~ ~ /xel^'^l ^ ~ /xel^^^l E,,^ ~ lO"" ^orb , (7) 

where i?orb is the energy (per unit mass) of the planet in the gravitational potential well of the 
star. As result, we estimate that Ep/Eo^]^ = ©(/ie'-^*'), so that the potential well of the resonance is 
about 10^ times shallower than the potential well of the star. As a result, planets can be removed 
from resonance much more easily than they can be removed from orbit (by turbulent fluctuations 
— note that planets are often ejected by scattering after orbit crossing). 



2.2. Turbulent Fluctuations 

The net effect of turbulence is to provide stochastic forcing perturbations on the mean motion 
resonance, which is modeled here as a pendulum. To include this stochastic process, the equation 
of motion ^ can thus be modified to take the form 

^ + [l + %5([T]-Ar)]sin0 = O, (8) 

where r]k sets the amplitude of the turbulent forcing and Ar sets the time interval. These quantities 
are derived in Paper I (including a discussion of where to introduce the turbulent term) and are 
discussed further below. 

Briefly, the time scale At is the time (in dimensionless units) required for the turbulence in 
the disk to produce an independent realization of the fluctuations. This time scale is typically one 
or two orbit times at the disk location in question (LSA04, Nelson 2005). Since the relevant part 
of the disk extends outside the orbit of the outer planet, and since the outer planet has twice the 
period of the inner planet, the time scale Ar ~ AijOq{2ti / Vt) = SvrtJo/^- 

Next we need to specify the forcing strength over the time interval Ar. In the original derivation 
of the pendulum equation (MD99) , one finds 

^ ^ ^ ^ f i^^l (9) 
df^ dt 3\JdtJ' ^ ' 

where the second approximate equality follows from the relationship between the mean motion and 
the orbital angular momentum. Converting to dimensionless form and integrating over one cycle 
to produce a discrete quantity, one finds 



(10) 
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The amplitudes [(AJ)/J]fe have been calculated for a variety of cases using MHD simulations (e.g., 
Nelson & Papaloizou 2003, 2004; LSA04; Nelson 2005). In basic terms, the torque exerted on a 
planet by the disk will be a fraction of the benchmark scale = 2TrGT,rmp, where S is the disk 
surface density (Johnson et al. 2006). The amplitude for angular momentum variations is thus 
given by A J = fxT iiT£,{Sti /tt), where the factor 87r/r2 is the time over which one realization of the 
turbulence acts, /t ~ 0.05 is the fraction of the benchmark scale Tjj that applies for a disk without 
a planet, and F/j ~ 0.1 is the reduction factor due to the production of a gap in the disk (and 
hence loss of disk material). Including all of these factors (sec Paper I), the relative fluctuation 
amplitude is given by [(AJ)/J]j.ms ~ 10"'^ under typical conditions. With this typical amplitude 
for [(AJ)/J]fc and Vl/ujq ~ 100 (see above), we expect that r^rms ~ 0.005. 

In general, the fluctuation amplitude will vary from case to case over a range of at least an 
order of magnitude due to different levels of turbulence and varying surface densities, or equivalently, 
varying disk masses. The amplitude quoted above is applicable for planets of Jovian mass. Lighter 
planets produce smaller gaps and have larger values of the factor and hence larger forcing 
amplitudes [(AJ)/J]/j. Note that even planets that are too small to clear a gap will still have an 
effective value of F/j < 1 because the wakes produced by the planet compromise the turbulent 
fluctuations in the vicinity (Nelson & Papaloizou 2004). In addition, different disk systems will 
retain their disk mass for a range of lifetimes, and this effect leads to another source of system 
to system variation. Finally, even for a given initial disk mass and disk lifetime, the amount of 
turbulence that resonant systems experience depends on how late (in the life of the disk) the planets 
are formed and/or captured into resonance (e.g., Thommes et al. 2008). These variations thus allow 
for a wide range of possible outcomes. 

For completeness, we note that turbulence is not the only possible source of stochastic fluctu- 
ations acting on planets in mean motion resonance. For example, a residual disk of planetesimals 
can provide a granular background gravitational potential and act in a qualitatively similar man- 
ner (Murray-Clay &l Chiang 2006). The formalism developed here can be used with any source of 
stochastic fluctuations with a given amplitude 77^ and forcing time interval Ar = u)Q{/S.t). 



2.3. Phase Space Distribution Functions 

An important technique used to analyze the behavior of stochastic differential equations — 
including those describing mean motion resonance — is to work in terms of phase space variables. 
For the stochastic pendulum considered here, we rewrite the equation of motion in the form 

^=V and ^ = -[l + r,fc<5([r]-Ar)]sin<^. (11) 

The probability distribution function P(r, V) for these phase space variables obeys the Fokker- 
Planck equation (e.g., Mallick &; Marcq 2004, hereafter MM04) which takes the form 

dP ,BP . ,5P 1 . 2 ,52p 
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As outlined in Paper I, the phase space diffusion constant D is specified by the ampUtude 77^ of 
the turbulent fluctuations and the time interval At required for them to attain an independent 
realization, i.e., 

ivl) 



D 



At 



(13) 



In most cases of interest, the angle (j) varies more rapidly than the velocity V (see Paper I and 
MM04). As a result, we can average the Fokker-Planck equation (jl2p over the angle <j) to obtain 
the corresponding time evolution equation for the (/)-averaged probability distribution function 
P{t,V). Throughout this paper, we use the same symbol (here, P) for the distribution function 
both before and after an averaging procedure; this choice simplifies the notation, but could result 
in some ambiguity, although the relevant version of the distribution function should be clear from 
the context (this issues also arises in §4). After averaging, the resulting Fokker-Planck equation 
([12]) becomes a basic diffusion equation 



dP 
97 



Dd'^P 



which can be solved to obtain the result 



P{t,V) 



{ttDtY/^ 



exp 



Dt 



(14) 



(15) 



This solutions shows that the energy of the pendulum can grow large at long times. In this limit, 
the kinetic energy dominates the potential energy term, and we approximate the energy using 
E ~ The resulting probability distribution function for the energy then can be written in 

the form 



PiE,T) 



ttDt 



1/2 



exp 



2E 
1>? 



(16) 



For this solution for the distribution function, the expectation value of the energy grows linearly 
with time in the long time limit, i.e., {E) = Dt/4. As the mean energy increases, the probability 
of the system remaining bound in a low energy (resonant) state decreases. Here we use > 1 as 
the requirement for the libration angle to circulate and hence for the resonance to be compromised 
(see Paper I for further discussion). The probability Pb of remaining bound is thus given by 



ttDt 



1/2 



dE 



exp 



2E 
15^ 



e-' dz = Erf(zo) 



(17) 



where zq = {2/Dt)^^'^ and where Erf(z) is the error function (e.g., Abramowitz & Stegun 1970). In 
the limit of late times, when z <^ 1, the error function has the asymptotic form Erf(z) ~ 2zl \pK, 
and the probability that the planetary system remains in resonance is given by 



^'b(r) 



1/2 



1/2 



where this expression is valid for sufficiently late times when Dt ^ 1. 



(18) 
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2.4. Numerical Integrations 

The third and final technique that we use to study this problem is direct numerical integration 
of the planetary systems. In other words, we integrate the full set of 18 phase space variables 
of the corresponding 3-body problem (using a B-S integration scheme, e.g., Press et al. 1992). 
Turbulence is included by applying discrete velocity perturbations at regular time intervals; for the 
sake of definiteness, the forcing intervals are chosen to be twice the orbital period of the outer planet 
(four times the period of the inner planet). Both components of velocity in the plane of the orbit 
are perturbed randomly, but the vertical component of velocity is not changed. The amplitude of 
the velocity perturbations is then allowed to vary. 

For the numerical experiments carried out for this paper, we start a large ensemble of A'ens 
planetary systems in the 2:1 mean motion resonance, and then integrate forward in time using the 
velocity perturbations (as described above) to model turbulent fluctuations. Throughout this work, 
we use the resonance angle defined by equation ([42]) for interactive systems such as GJ876 (see §4), 
and its analog in the circular, restricted 3-body problem given by equation ([T]) with {j^l = 2. Our 
numerical experiments show that this choice of resonance angle allows the systems to remain in 
bound (resonant) states for longer times than for the standard case where |j4| = 1 (see Paper I and 
§4). Since we are primarily interested in how turbulent fluctuations can compromise mean motion 
resonance, we want to consider the cases that are most robust (longest-lived). Finally, we note that 
the analysis of this paper can be performed for a host of other forms for the resonance angle. 

Over the course of time, systems move both in and out of resonance. To monitor this behavior, 
we determine the maximum libration amplitude of the resonance angle over a fixed time interval. 
This time scale is chosen to be somewhat longer than the libration time, which is much shorter 
than the evolution time, i.e., the time scale over which a substantial fraction of the systems leave 
resonance. Note that the libration time is typically ~ 100 orbits (see equation [3]), whereas the 
evolution time is of order 10^ orbits. For the sake of definiteness, we consider a planetary system 
to be in resonance when its maximum libration angle is less that 90 degrees; for larger maximum 
libration angles, we consider the resonance to be compromised. Note that the boundary at 90 
degrees is somewhat arbitrary; fortunately, however, the statistics for the fractions of bound systems 
is insensitive to this value (see Paper I). 

Although the observed set of extrasolar planetary systems shows an astonishing degree of 
diversity, here we consider only two particular examples that bracket the possibilities. For the 
first case, we consider analogs of the GJ876 planetary system, which has two planets that are 
observed to be deep in a 2:1 mean motion resonance (Marcy et al. 2001). This system has a 
relatively small central star (with mass of only M* = 0.32 Mq) and relatively large planets (with 
masses m = 0.79 and m2 = 2.52 mj). The eccentricities are substantial, with e = 0.26 for the 
inner planet and 62 = 0.034 for the outer planet. The periods are approximately 30 and 60 days, 
respectively (see Rivera et al. 2005 for further detail). With these parameters, the GJ876 system 
is the most highly interactive system observed to date. Although comparable solar systems with 
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larger planetary masses and/or higher orbital eccentricities would be even more interactive, they 
would also most likely be unstable (Laughlin & Chambers 2001). Since the pendulum model for 
mean motion resonance represents the simplest possible model, and in particular does not allow for 
interactions between the planets, the GJ876 system (its analogs considered here) is about as far as 
possible from the idealized one variable pendulum model. Because of its highly interactive nature, 
the GJ876 system analogs take a long time to integrate and we consider ensembles of Mens = 10^- 

At the other end of parameter space, we consider a model solar system with relatively little 
interaction between the two planets. In this case, we take the outer planet to be a Jupiter mass 
planet (m2 = Imj) in a 150 day orbit around a 1.0 Mq star. The second, inner planet is taken to 
be a "Super-Earth" with mass m = O.Olmj. The period of the inner planet is taken to be 75 days, 
and the system is started in the 2:1 mean motion resonance. The orbital eccentricities are e = 0.1 
for the inner planet and 62 = 0.01 for the outer planet. This latter system is thus well approximated 
by the circular restricted 3-body problem, which is used in the derivation of the pendulum model 
for mean motion resonance (DM99). We thus expect ensembles of this class of solar system to 
follow the predictions of the stochastic pendulum model. For these systems, the integrations can 
be run faster and we consider ensembles of A^ens = 10^. 

The effects of turbulence on mean motion resonance are illustrated by Figures 1-3. The first 

two of these figures show the fraction of resonant states as a function of time for an ensemble of 
solar systems that started out with the configuration of GJ876, and Figure 3 describes results for 
the less interactive Super-Earth systems. 

Figure 1 shows the various ways to account for the fraction of bound states. The solid curve 
shows the ratio of the number of planets in resonance to the number of planets in the original 
sample. However, some planets are lost with time: systems that leave resonance often experience 
orbit crossings, which can lead to planetary ejection. The dotted curve shows the fraction of systems 
that have not ejected a planet. Finally, the dashed curve shows the ratio of the number of planets 
in resonance to the number of planets that remain in orbit. For the remainder of this paper, 
we use this latter ratio (number of systems in resonance over the number of surviving systems) 
as the relevant fraction of bound states (this ratio is the most directly observable). Note that 
a substantial fraction of the systems that leave resonance lose a planet. This ejection takes place 
because of orbit crossing and subsequent planet-planet scattering. Systems are protected from such 
behavior while they remain in resonance, but interactive systems (e.g., GJ876) readily lose planets 
in when resonance is compromised. 

Figure 2 shows how the time evolution of the bound fraction depends on the level of turbulence. 

For an ensemble of GJ876 systems, the figure shows the fraction of bound states as a function of 
time for two levels of turbulent fluctuations, where {Av)/v = 0.0002 and 0.0004. Note that in 
both cases, the bound fraction is nearly a straight line on this log-linear plot, so that the time 
evolution is nearly exponential with a well-defined decay rate, or, equivalently, half-life (see §4 for 
further discussion). Notice also that the diffusion constant D sets the time scale for systems to 
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Fig. 1. — Results from numerical integrations of analogs of the GJ876 system to illustrate varying 
definitions of the bound fraction. Here the turbulent fluctuations have amplitude {Av)/v = 0.0002. 
The solid curve shows the ratio of the number of planets in resonance to the number of planets in 
the original sample; the dashed curve shows the ratio of the number of planets in resonance to the 
number of planets that remain in orbit; the dotted curve shows the fraction of planets that remain 
in orbit. 
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Fig. 2. — Time evolution of the bound fraction (number of planets in resonance divided by the 
number of planets that remain in orbit) for numerical integrations of the GJ876 system. The two 
solid curves show the numerical results obtained using two levels of turbulent fluctuations, with 
{Av)/v = 0.0002 (right) and 0.0004 (left). Note that the diffusion constant D sets the time scale 
for systems to leave resonance, and that D depends on the square of the fluctuation amplitude, so 
the decay times for the two numerical experiments should differ by a factor of four. For reference, 
the two dashed lines show purely exponential decay with time scales that differ by a factor of four. 
The error bars show the expected amplitude of root-AT fluctuations, which are somewhat smaller 
than the observed variations in the numerical results. 
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Fig. 3. — Time evolution of the bound fraction for numerical integrations of the Super-Earth 
system, where m2 = nij, P2 = 150 days, m = 0.01 mj, and = I.QMq. The bound fraction 
is defined to be the number of planets in resonance divided by the number of planets that remain 
in orbit. The two solid curves show the numerical results obtained using two levels of turbulent 
fluctuations, with {Av)/v = 0.0001 (top) and 0.0002 (bottom). At late times, when the bound 
fraction falls below P^, < 0.1, the dashed lines show the power-law behavior Py, t^^/^ predicted 
by the analytic model (which is valid in the absence of planet-planet interactions). The error bars 
show the expected amplitude of root-A^ fluctuations. 
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leave resonance, and that D depends on the square of the fluctuation amphtude, and hence the 
decay time scales differ by a factor of four. The dashed lines shown in Figure 2 show a purely 
exponential decay, with time scales that differ by exactly a factor of four. The agreement shown 
in Figure 2 thus confirms our theoretical expectations. The error bars plotted on the dashed lines 
show the expected amplitudes of root-A^ fluctuations, which are somewhat smaller than the actual 
variations found in the numerical results; the presence of planet-planet interactions (see also §4) 
increases the level of these variations. 

Figure 3 shows the time evolution of the fraction of systems remaining in resonance for an 
ensemble of Super-Earth systems. In this set of experiments, two different levels of turbulence are 
used, so that the amplitude of the velocity perturbations are {Av)/v = 0.0001 and 0.0002. This 
class of systems is close to the idealized, circular restricted 3-body problem, and hence the behavior 
of the resonance is expected to follow the predictions of the stochastic pendulum model. As shown 
in Figure 3 (by dashed lines) the fraction of bound states decreases as a power-law in time, with 
Pb ^ t~^^^ at late times, roughly as predicted by §2.3. Given that the velocity perturbations 
differ by a factor of two, one expects the coefficients of the power-laws to also differ by a factor 
of two. However, the numerical results show a somewhat wider gap, namely a factor of ~2.6. 
Although the size of the leading coefficient is correct in order of magnitude, the following three 
ambiguities arise in determining its predicted value from the analytic formulation: The first issue 
is the mismatch between the simplified treatment of the circular, restricted 3-body problem and 
full (18 variable) integrations; as a result, the libration frequency of the numerical system is not 
exactly the same as that given by the approximation of equation ^ . The second difference is due 
to the presence of a partial barrier; after systems leave resonance, they have a somewhat lower 
probability pret of re-entering resonance from an excited state because the physical system (with 
18 phase space variables) is more complicated than the one-variable pendulum (see Paper I for 
further discussion). This effect should not be overly severe for these Super-Earth systems, as they 
are close to the regime of validity of the pendulum approximations (MD99), but the partial barrier 
is still present. Yet another relevant process is that planetary ejection can take place. In 3-body 
numerical integrations of Earth-like planets with Jovian companions, when the periastron is the 
same as in these systems, the expectation value of the ejection time is about 3 x 10^ orbits, or 
6 X 10^ yr (David et al. 2003). In spite of these complications, however, the basic trends shown by 
the numerical results are in qualitative (and rough quantitative) agreement with the expectations 
of the simplest pendulum model. Finally, we note that the error bars shown in Figure 3 show the 
amplitude of root-A^ fluctuations, which are roughly the same amplitude as the variations seen in 
the numerical results. 

Now we can compare Figure 2 with Figure 3: For the interactive GJ876 systems of Figure 2, the 
fraction of bound states Pb is presented as a log-linear plot, so that straight lines (as found in the 
simulations) correspond to exponential decay. For the Super-Earth systems of Figure 3, however, 
the fraction Pb is presented as a log-log plot, where straight lines (as found in these simulations) 
correspond to power-law decay. Since exponential decay is much stronger than power-law decay. 
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we expect that highly interactive systems will leave mean motion resonance much more readily 
than less interactive systems. This effect is accounted for in the analytic model developed §4. In 
particular, this model predicts power-law decay of the fraction of bound states in the regime of 
minimal interactions between planets, and exponential decay in the regime of strong interactions. 



3. STOCHASTIC PENDULUM WITH DAMPING 

This section considers the effects of damping on the maintenance of mean motion resonance 
in the face of turbulent fluctuations. We begin with the dimensionless equation for a stochastic 
pendulum (see equation [6]), include the turbulent forcing term (see equation [8] and Paper I), and 
then add a damping term, 

^ + 7^ + [1 + VkS{[r] - At)] sin = , (19) 

where 7 is the damping rate for the resonance angle. Here we take 7 to be constant; its value 
is not well determined, but we expect it to be small. For example, during planet migration, the 
semimajor axes decrease with an effective "damping rate" 7a = —d/a and the disk tends to damp 
the eccentricity at a rate je = — e/e. Since planet-planet interactions excite eccentricity during the 
migration epoch, the net damping rate 7 of the resonance is generally much smaller than either 
7e or 7(j. As one example, even in the extreme case of eccentricity damping where je = 1007a, as 
required to explain the observed orbital elements of GJ876, the libration amplitude of the resonance 
angle remains nearly constant, so that 7 ~ (see Figure 4 of Lee & Peale 2002). 

The basic equation of motion ()19p can be converted into a system of first order differential 
equations by introducing V, i.e., 

^ = V, ^ = -jV-[l + %<5([r] - At)] sin . (20) 
ar CLT 

This set of equation corresponds to the following Fokker-Planck equation 

dP /„ ..dP\ . dP D . 2 ,d^P 



where the effective diffusion constant D = {rjf.) /At determines the strength of the stochastic forcing. 
The terms that contain 7 encapsulate the effects of damping. As the next step, we use the fact that 
(j) changes more rapidly than V , and hence we can average out the cf) dependence (see §2, Paper I, 
and MM04). This averaging procedure removes the sin^ and d/d<j) terms, and the Fokker-Planck 
equation simplifies to the form 

dP Dd'^P ^BP 

This equation can be considered as a "modified" diffusion equation. Specifically, we can convert 
this equation into a diffusion equation in standard form by making a suitable change of variables. 
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I.e., 



T = 



D 

87 



1)> 



U = Ve^^ 



P = Qe 



With these definitions, we note that 



dP 



yrdQ 

dr 



dQ dT dQ dU 



dT 



dU 



and that 



dP 



dQ 



and 



d^p 



d'^Q 



dV ^ dU 9[/2 • 

Using these results in the original equation, we find a diffusion equation of the form 



P7 + e-,C/j^+^_ 
which then reduces to the simpler form 



dt 



dm 



rdQ 

dU 



dQ ^ d^Q 
dT ~ ■ 



(23) 

(24) 
(25) 

(26) 
(27) 



This latter result is thus an ordinary diffusion equation. The initial condition considered here is 
given by P{t = 0;V) = S{V), so that all of the oscillators start in a bound (resonant) state. In 
terms of the new variables, this condition corresponds to Q{T = 0;U) = S{U), and the solution for 
the distribution Q{T, U) can be written in the form 



1/2 



87 



1/2 



exp 



D(e27^ - 1 

Finally, the distribution in terms of the original variables is given by 



27y^e 



2„27r 



D (e27^ - 1) 



(28) 



P{t,V) 



2je 



2jT 



L»7r(e27^ - 1) 



1/2 



exp 



27y^e 



D (e27T _ 1) 



(29) 



Given the analytic form of equation (|29p . we can now consider the asymptotic forms. For the 
limit where 7r ^ 1, for small damping and/or short times, we can expand the exponentials to first 
order and find 

Pi^^ V) = / ^ M/o exp 



I!' 

Dt 



(^Z)r)V2--- ' 

which is the same result obtained earlier (see equation [15] and Paper I) with no damping. In the 
opposite limit where 7r ^ 1, the distribution takes the form 



P{t,V) 



_27 
Dtt 



1/2 



exp 



271/2 
D 



(31) 



i.e., the distribution approaches a gaussian form that is constant in time. 
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The general behavior described here can be understood through the fohowing heuristic argu- 
ment: In the equation of motion (j20p for V, the time derivative has two contributions. For the 
stochastic part, V"^ ~ Dt so that V ~ VDt and hence V ~ V/2t. For the damping contribution, 
V ~ —"jV. The two terms have opposite signs and are equal in magnitude when 2jt = 1. For 
smaller values of 7T, the stochastic term dominates and we recover the results of the undamped 
stochastic oscillator. For larger values of jt, the damping term balances the stochastic term and 
the distribution P{V) becomes stationary. 

Given the distributions derived above, the expectation value of E = the kinetic energy 

of the oscillators, is given by 



1 



D e^T^ - 1 



87 



D 

87 



[1 



-27Tl 



(32) 



In the limit of no (small) damping, the energy expectation value becomes {E) Dt/4:, in agreement 
with the result from the undamped case (Paper I). However, for any nonzero value of 7, the 
expectation value approaches a constant {E) D/(8'j) in the long time limit r 00. 

Numerical simulations of the stochastic pendulum, analogous to those presented in Paper I, 
illustrate the behavior indicated by equation (j32p . Specifically, Figure 4 shows the energy expec- 
tation value for an ensemble of stochastic pendulums as a function of time. The two solid curves 
show the result of numerically integrating the stochastic pendulum equations (see also Paper I) for 
fixed diffusion constant D and two values of the damping rate 7 (here, in dimensionless units, D = 
2 and 7 = 0.005, 0.05). The two dashed curves show the prediction of equation ([32]) for the same 
fluctuation amplitude and damping rates. The two treatments are in excellent agreement. At early 
times, the energy expectation value increases linearly with time, but then saturates and approaches 
an asymptotic value (-E)oo = D/{8j) at late times. 

Next we consider the bound fraction, i.e., the fraction of the systems that remain in resonance 
as a function of time. For the sake of definiteness, we define the bound fraction to be the 
fraction of states with E = < 1, i.e., we consider only the kinetic energy in the accounting. 

The bound fraction is thus given by the integral 



After defining new variables 



27 



o7T 



exp 



2-iV^e 



2„27r 



L>(e27T _ 1) 



dV. 



2-iV^e 



2p27r 



D (e27T _ 1) 

this integral can be simplified to the form 



and 



476^^^ 



D (e27r _ 1) ' 



^0 



71" Jo 



Erf(zo) 



(33) 



(34) 



(35) 
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Fig. 4. — Energy expectation value for an ensemble of oscillators as a function of time. The two 
solid curves show the result of numerically integrating the stochastic pendulum equations for fixed 
diffusion constant (D = 2 in dimensionless units) and two values of the damping rate (7 = 0.005 
and 0.05). The two dashed curves show the predictions of equation (|32p . 
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where Erf(2;) is the error function (e.g., Abramowitz & Stegun 1970). We are often interested in 
the limit where the argument of the error function is small, so that we can expand Erf(z) in a 
power series in z: 



Keeping only the first order term, we can write the fraction of bound states in the form 

/ ^ \l/2 pIT- 

This expression is valid in the hmit ^ li which requires the following conditions to be met 

» 1 and L* > 47 . (38) 

The first condition will be satisfied whenever turbulent fluctuations persist for a long time, as 
expected in circumstellar disks associated with young stars; this requirement is necessary for the 
bound fraction to be small in the absence of damping (Paper I). The second condition requires 
the diffusion constant to be larger than the damping rate. In other words, the size of the ratio D/7 
determines whether turbulence or damping dominates the dynamics. 

Now we can consider the limiting forms of the result from equation (|37p . In the limit of small 
damping or sufficiently short time scales, 7r <C 1, we recover the prediction for the bound fraction 
from the undamped calculation of Paper I, i.e., 

limPb=^ . (39) 
7T^o \ttDt J 

In the opposite limit where 7r 3> 1, the bound fraction approaches an asymptotic (constant) value: 

hm P^ = A(^y^^ . (40) 

The asymptotic value of Pb will thus be small provided that the damping rate 7 is small compared 
to the diffusion constant (see equation [38]). 

Figure 5 shows the results from integrations of two ensembles of stochastic oscillators; the 
resulting behavior is in good agreement with the analytic predictions derived above. This Figure 
shows the fraction of bound resonant states as a function of time for the same two ensembles of 
systems considered in Figure 4, i.e., stochastic pendulums with turbulent fluctuations and two 
different values of the damping rate 7. As expected, the bound fraction initially decreases like a 
power-law as indicated by in equation ([39]) . where Ph{t) ~ and eventually saturates at the 

value given by equation ([^0]) . 



In order to estimate the size of this asymptotic value for Py,, we rewrite the limiting expression 
ni) in the form 

'?rms (A^orb) ^^^^l^orb) ' ^^^^ 
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Fig. 5. — Fraction of bound (resonant) states for an ensemble of stochastic oscillators as a function 
of time. The two solid curves show the result of numerically integrating the stochastic pendulum 
equations for fixed diffusion constant (D = 2 in dimensionless units) and two values of the damping 
rate (7 = 0.005 and 0.05). The two dashed curves show the predictions of equation (j37p . The 
error bars plotted at r = 700 show the expected amplitudes of root-A^ fluctuations (which roughly 
bracket the variations found in the numerical results). 
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where N-y = is the number of damping times and A'orb is the number of orbits for which 
turbulence is active (see also the definition of equation [13] ) . Unfortunately, the relevant value of the 
damping time scale remains uncertain. The time scale for planetary migration - the time required 
for the semimajor axis a to change - is typically a few Myr in these systems, roughly comparable 
to the disk lifetimes. In addition to driving migration, the disk tends to damp eccentricity (for the 
most part - see Moorhead & Adams 2008, Goldreich & Sari 2003) on a time scale that is roughly 
comparable to the migration time (e.g., Crida et al. 2008). If we assume that turbulence is active 
over the same time interval when damping is active, a few Myr, then ~ 3 and the number of 
orbits A^orb ~ 10^. These values would predict an asymptotic value of the bound fraction Pb ~ 0.2. 
However, although eccentricity damping leads to a damping of the resonance excitations, planet- 
planet interactions lead to eccentricity excitation and tend to counteract its effects. As a result, 
values of 7 and are smaller than in this simple estimate, and the asymptotic value of Pb will 
also be smaller. Further, all of the relevant parameters (the size of the damping, the duty cycle 
of the turbulence, the amplitude of the fluctuations, etc.) are expected to vary from system to 
system. We thus expect damping effects to counteract turbulence in some cases, but not in others. 
In any case, the results of this section determine how the asymptotic value of bound fraction Pb 
depends on the amplitude and duty cycle of the turbulence, and the damping rate 7. 



4. RESONANCE MODEL INCLUDING INTERACTIONS 

The analysis carried out thus far indicates that the fraction of resonant states tends to decrease 
exponentially with time for highly interactive systems (e.g., GJ876) in the presence of turbulence. 
On the other hand, the fraction of bound states in less interactive systems (e.g., the Super-Earth 
system) decreases as a power-law in time, in agreement with the pendulum model of resonances 
(with stochastic fluctuations). The pendulum model, with only one variable, allows systems to 
freely random walk in and out of resonance; highly interactive systems seem to have more difficulty 
re-entering a resonant state. These results suggest that the barrier to returning into resonance is 
due, in part, to interactions between the two planets. These same interactions also lead to such 
systems leaving a resonant state more easily in the presence of stochastic fluctuations. 

In this section, we derive an equation for the evolution of the resonance angle where interactions 
between the planets are included (see also Holman &; Murray 1996, Quillen 2006). As shown below, 
this analysis initially retains many variables, and thus includes more physics than the (one- variable) 
stochastic pendulum considered previously. In order to obtain tractable results, however, we average 
over many of the variables and find a Fokker-Planck equation for the reduced and averaged problem. 
For the solution to this Fokker-Planck equation, the fraction of bound (resonant) states decreases 
exponentially with time when planet interactions are important; when interaction terms are small, 
the solution leads to the same power-law behavior oc r"^/^ found previously (see equation |18j 
and Paper I). 
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4.1. Equations of Motion 



In this section we derive the equations of motion to describe the evolution of one particu- 
lar resonance angle and the corresponding orbital elements. Note that we cannot use the same 
resonance angle as in §2. In that case, the simple pendulum equation resulted from the circular 
restricted three-body problem (see MD99). In this context, we want to allow for interactions be- 
tween the planets (so we must move beyond the restricted three-body problem) and hence the outer 
planet can attain nonzero and time varying eccentricity (so the circular approximation is no longer 
applicable). We thus consider the following (non-standard) resonance angle 



2A2 — A -|- 2tu2 — 2uj , 



(42) 



where A is the mean longitude and w is the argument of periastron. Note that this resonance angle 
is the analog (beyond the circular, restricted 3-body approximation) of that from equation ([1]) for 
the case |j4| = 2. Throughout this derivation, the variables associated with the outer planet are 
denoted with the subscript '2', whereas the variables associated with the inner planet are left with 
no subscript. 

For our chosen resonance angle (I42p . the lowest order part of the disturbing function for the 
inner planet is given by 



(7^) 



Gm2 
02 



7^(^sec)^^2g2 [/,(„) + /^(a)] cos, 



and that for the outer planet is 



(7^2 



Gm 



^^(sec) ^ ^2g2 \^afd{a) + fi{a)] cos ( 



(43) 



(44) 



where e is eccentricity, m the mass, (0,02) are the semimajor axes, and G is the gravitational 
constant. The functions fd-, fe, fi are defined and discussed in MD99. To lowest order, the direct 
secular contribution is then given by 



7^ 



(sec) 



D 



(e^ + ei)/s,i(a) + 662/^,2(0) cos(tU2 - w). 



(45) 



The corresponding equations of motion for the orbital elements of the inner planet then become 

n = — -— = ^ e^ei [fd{a) + /e(a)] smc/), (46) 



a2 dX 
1 dn Gm2 



0^02 



[ee2/s,2(a) sin(tU2 - w) + 2e^e| [/d(a) + /e(a)] sin , (47) 



■uj 



1 dTZ- G'm2 

^^2g gg = na?a2e [^^-^^'i*^") + 62/5,2(0) cos(ti72 - w) + 2eel [fd{a) + fe{a)] cos cf] . (48) 

To complete the set, note that e = e'^'Uj/2, where e is the mean longitude at epoch. However, 
because e is times smaller than w, we can ignore time derivatives of e in favor of time derivatives 
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of zu. This assumption is valid within the range of validity for these equations, which are the lowest 
order terms in an expansion in the eccentricities (6,62). The equations of motion for the orbital 
elements of the outer planet are similar, i.e., 



^ 3 dTZ2 6Gm 9 2 r ^ / \ , ^ / m ■ / /^r,\ 
n2 = 2 = -^—^ ^2 + smcj), (49) 



^2 = ^ [ee2afs,2ia) sm{w2 - w) + 2e'^el [afaia) + fi{a)] sin cp] , (50) 

\l2a2ae2 

^2 = [2e2Q;/s,i(Q:) + eafs,2{oi) cos(ro2 - w) + 2e^e2 [afd{a) + fi{a)] cos ^1 , (51) 

Next we rewrite X as \ = Clt + e, and assume that tit <^ $7, and that e ^ zu. The time 
evolution of the resonance angle is then given by 

4)^2^2-^ + 2tU2 - 2w . (52) 

Also, one can write GM* = Sl^o^ = i^lal, where M* is the mass of the central star. Using this 
latter result, we define 

Cr ^ §^ [h{a) + /.(a)] = ^Sla [/j(a) + h{a)\ . (53) 

(2) ^21 (2^ 

and define Q , ™ analogous way (where the superscript now refers to the second 

planet). With this choice of notation, the equations of motion take the form 

J7 = -SC^Oe^e^sin^, (56) 

e = — Cs2e2 sin 6 — 2Cree\ sin 0, (57) 
62 9 

W = 2Csl+Cs2 — cos + Cre2 COS ^, (58) 

$72 = 6C^J72e^eisin(^, (59) 
62 = Cg^e sin 6 + 2C^^^e^e2 sin (60) 
tt72 = 2C2^ + C^^ — COS e + e2 cos 0, (61) 

62 

where 6 = W2 — w. 

After taking a second time derivative of w and W2^ we can derive the equation of motion, 
which gives (p in terms of the variables (^, e, 62, fi, ri2) ^)- For simplicity, we also leave the variables 
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(e, 62, 0) in the equation, but these can ah be rewritten in terms of the variables Usted above. The 
resulting equation of motion then takes the form 



36^62 f44^)j72 + CrQ) sin(/. + 2 



^(2) e 

62 



1 + 



Cg262 



+ C,2- 1 + 



Cs262 



COS ( 



+ 2 C 



^(2) 



•'s2 



62 



s2- 



Cr62e2 ) COS</) + 2 [C. 



'{2)g2 



)sm( 



(62) 



^sin0 + 2 (Cp)66 

62 6 

Consider the ordering of the terms in equation (j62p . All of the time derivative terms (i.e., those 
including e or </>), have an extra factor of C, which introduces another factor of /U = m/M*. All 
of the terms in the equation are thus second order in m /M* , except those coming from the Vt and 
^2 terms, which are only first order. As a result, to leading order, we recover the usual pendulum 
equation of motion to describe a mean motion resonance, i.e., 



o 2 2 
36 62 



^^^'^^2 + Cr^ sin 



(63) 



where the effective frequency of the system is a linear combination of the frequencies of the "oscilla- 
tors" for each planet separately. Notice that the frequency of this pendulum contains extra factors 
of eccentricity compared to the simplified case considered in §2; this difference arises because the 
two pendulum equations ^ and (|63p correspond to different resonance angles. 



4.2. Inclusion of Turbulence 

Now that we have derived the "classic" equation of motion (j), without fluctuations, we must 
include the effects of turbulence. In this treatment, we incorporate turbulent forcing as a set of 
discrete impulses in the equation of motion for the eccentricity, i.e., 

6 = — Cs262 sin0 — 2Cr6e2 sin0 + ^ where ^ = rjk5 {[t\ — At) . (64) 

As shown above, the equations of motion for all of the variables include eccentricity, and hence 
this ansatz for including turbulence is convenient. These impulses acting on the eccentricity then 
produce corresponding impulses in the equation of motion for the resonance angle (p, and in the 
equations for zu and zu2- Since the orbital elements of the external planet should not be directly 
perturbed by an impulse acting on the inner planet, the validity of including impulses in the latter 
equation is somewhat subtle. However, the effect shows up in the second derivative, which acts 
like a force, which in turn changes when the inner planet experiences a perturbation (here due 
to turbulence). Finally, we note that the distribution of impulses produced by equation ([64|) is 
more complicated than the uniform distribution of velocity perturbations used in our numerical 
simulations. Because of this complication, the amplitude of the eccentricity perturbations will be 
different from the effective amplitude of the perturbations acting on the other variables (e.g., the 
velocity V = (j)). 
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We can now write out the Fokker-Planck equation, which includes six variables at this stage 
of derivation. In order to obtain tractable results, we average over most of the oscillating variables. 
To start, we note that the equation simplifies considerably after we average over the angle (p 
itself. Although this approximation is often used (e.g., MM04), its implementation in this case is 
somewhat complicated. Here, the frequency of the largest oscillations of (f) (those lowest order in 
C) is proportional to VC, whereas the primary frequencies of the other quantities are proportional 
to C; note that C <C 1, so that (p oscillates with higher frequency. Further, the terms for which this 
frequency difference does not hold are small and can be neglected. As a result, (/) has effectively a 
higher frequency and can be averaged out (sec Paper I and MM04 for further justification). After 
averaging, the Fokker-Planck equation becomes 



62 ^ 



62 e 



'sin(ti72 — vj)—— + 



dP v. d^p . ^(2) 1 f ^ , c,2el \ . d'p 



dV 2 de^ 



62 



1 + 



^S2 ^ y 



cosO 



dVde 



Ci?- 1 + %^ -s^ 



d^p 



(65) 



Because of the difference in perturbation amplitudes for the eccentricity e and for the velocity 
V, the diffusion constants Vg, V^y, and are not equal (in general). 



4.3. Apsidal Resonance 



Next we consider how the angle 9 evolves with time. We start with the equation of motion for 
9 in the form 



e 
62 



I^s2 + -2 



cos^+ ( Cfj—-Cs2— ] 9sm9 

62 ^ 



+ 2 (^C^^^ee + Cre2e2) cos ^ + 2 (cl^'^h'^ - Crei) ^ sin ( 



(66) 



As discussed above, in this approximation we average over the (rapidly oscillating) variable (p. This 
averaging reduces the equation of motion for 9 to the form 



9 



62 



- (ci')+c,2^ \ --iCs2+cii;'^ 



's2 "^'-s2-2 



62 



^(2) 



=2\ H 



's2 -r <^g2 -f 



cos ^ + ( Cf^ — -Cs2—]9sm9, (67) 

62 ^ 



where now e, 62, ^ are ^-averaged, and thus obey the reduced equations of motion 

6 = -Cs2e2sin6', (68) 
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e2=c2^esin0, (69) 
(2) o^, , /^^(2)e ^ ^2 



= 2C^^^ - 2Csi + (^C^2 - - Cs2^ J cos 0. (70) 

When the masses of the planets are comparable, then C ~ C^'^\ and the first two terms on the right 
hand side of equation ()70p nearly cancel. If, in addition, the planets have the appropriate values 
of eccentricity (e, 62), then the final terms can also nearly cancel out. The small remaining part of 
the first two terms can then balance that of the second two terms, ~ 0, and the system can be 
trapped in a "^-resonance" or apsidal resonance (for further discussion, e.g., see Chiang & Murray 
2002, Murray-Clay & Chiang 2006). Numerical experiments show that the Super-Earth systems 
considered in this paper are not in apsidal resonance, whereas the analogs of the GJ876 system are. 
Further, if the planetary masses in the GJ876 systems are changed slightly, then the value of e/e2 
adjusts so that « 0. 

To see how this latter adjustment works, consider the time derivative of e/e2, i.e., 

-- — = 2" = -Cs2sm6' - cy-g sm6'. (71) 

dt \e2j 62 6^ 

(2) 

Keep in mind that both Cs2 and C).2 are negative. Now, suppose that 6 is increasing with sin0 > 0. 
Then e/e2 would increase. This behavior makes the positive term — Cs2(62/e) cos 9 in equation (|70|) 

(2\ 

decrease in magnitude, and the negative term, C),2 {e/e2) cos 6, increase in magnitude, which would 
cause 9 to begin to decrease. The same argument can be used to show that e/e2 would decrease if 
9 became negative, which would result in 9 becoming more positive. 

When the system is not in apsidal resonance, the variables (^, e, 62) all oscillate at a frequency 
given by a characteristic value of 6. As shown by equation (j70p . this characteristic value of 6 is 
of the order C. Likewise, when the system is in apsidal resonance, so that ^ ~ 0, the equation of 
motion for 9 reduces to the approximate form 



r ,r ^2\ , ^(2) /^^ ^(2)6^^^ 
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(72) 

'2/ J 

In this case, the variables 9, e, and 62 still oscillate at a characteristic frequency given by C. This 
result justifies averaging over the angle (j), which has a characteristic oscillation frequency of \/C 
(keep in mind that \/C ^ C because C <C 1). However, we cannot average the Fokker-Planck 
equation over 9 (i.e., over w and W2) in either case, because e and 62 both vary on the same time 
scale and because the amplitude of this variation is substantial. 

Next we consider the effect of turbulence on 9. As above, we will add a turbulent forcing term 

to e, 

6 = — Cs2e2 sin6' + ^ (73) 

where ^ has the same meaning as before. We can now write a system of first order differential 
equations for 9,e,e2, 

6 = to, (74) 
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and find their equivalent Fokker-Planck equation: 
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Here we use Q to denote the distribution function for the apsidal resonance angle 6^ which is not to 
be confused with the distribution function P for the original resonance angle 0. This equation (j78p 
has almost the same diffusion terms (those containing the diffusion constants 2?e; ^ei^j and T)^^) as 
the Fokker-Planck equation (j65|) for P (only the factors of 2 are different). This correspondence 
suggests that for systems starting in apsidal resonance, the bound fraction for apsidal resonance 
should show qualitatively the same time dependence as that of the (/)-resonance. For example, we 
observe exponential decays in the bound fraction for both angles in the GJ876 system (see also 
Paper I). Furthermore, the effective diffusion constant for apsidal resonance is approximately four 
times larger than the diffusion constant for removing systems from (/)-resonance. Systems thus tend 
to leave apsidal resonance four times faster, which is a trend seen in our numerical simulations. 



4.4. Averaging the Fokker-Planck Equation 

To make further process, we need to average over 6. This must be done carefully because all 
the terms in the problem vary on the same timescale as 9. We assume that the system is not in 
^-resonance, and that 9 = ujt, for some w. As a result, both e and 62 can vary with substantial 
amplitude about their mean values (denoted here as {e) and (6)2) in the following way: 

e = (e) + ^^^cos0= (6) + <56cos6', (79) 

UJ 

e2 = (e)2- ^ ^ cosg= (6)2-^2 cosg, (80) 
u> 

where the second two equalities define 6e and 5e2- This ansatz assumes that these perturbations 
of the eccentricity are smaller than the eccentricities e and 62- However, these perturbations are 
nonetheless significant and, as derived below, provide the effects of interactions in this formulation. 
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From the original equations of motion, we see that to ~ 2C^^ — 2Csi- We can now average over 
9 in equation (j65p . and derive the simphfied form for the Fokker-Planck equation 
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where only the first order terms in 5e and 5e2 have been retained. In general, we can expand these 
formulae in a power series in Se/{e) and keep higher order terms. For simplicity, we use only the 
leading order terms here. When 62 <C e, we can use the result 
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and we can evaluate the integral (over 9) exactly by contour integration. In this case, we find 
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As a consistency check, note that if we reduce equation (|81|) using the approximation of equation 
(f82]) and then evaluate equation ([83]) in the limit 662/ {e) 2 ^ 1, the resulting expressions agree. 

At this point, it is advantageous to average over the eccentricity e. In the Fokker-Planck 
equation, the second derivative terms with respect to e suggest that the eccentricity is undergoing 
an ordinary diffusion process, in which case its expectation value would grow like \/t. In the 
true physical system, however, the eccentricity is bounded both from below (e > 0) and from 
above (e < 1). The effective upper bound is even smaller because planets that reach a sufficiently 
large eccentricity often experience orbit crossings and eventual ejection. For sufficiently long time 
scales, the eccentricity is expected to approach a nearly stationary distribution, and it becomes a 
reasonable approximation to average over e. We thus define the following quantities: 
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where p{e) is the distribution of e. This approximation is straightforward when the original distribu- 
tion P{t, 6, V) is a separable function; when this is not the case, the meaning of the approximation 
is more complicated. In either case, after averaging, the Fokker-Planck equation becomes 
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which is now a one-dimensional diffusion equation that can be solved analytically using standard 
methods. 



4.5. Solution and Results 



Before proceeding further, we simplify notation by defining constants {X, Y, Z) so that the 
averaged version of the Fokker-Planck equation (derived above) can be written in the form 



f = ^ + >'l^ 

dt dV 
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Using standard methods, we can write the corresponding solution in the form 
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This solution corresponds to a distribution which is diffusing with diffusion constant Z, drifting like 
Yt, and decaying on the timescale 1/\X\ (note that X is expected to be negative to leading order). 
In practice, however, we find that the quantity X is small; our numerical simulations indicate that 
V^/'Df, ~ Iff^ in typical cases. This ratio of diffusion constants is the same as the ratio of the 
potential energy of the resonance to the potential energy of the planetary orbit (equation [7]). The 
diffusion constant measures the efficacy of turbulence in changing the eccentricity, i.e., changing 
the orbit; the constant 2?^ determines how easily turbulence can change the speed V = (p, which 
is easier to change by a factor of ~ 10^. Because 2?e is small, the decay term in equation (j88p can 
be ignored on the intermediate time scales relevant to this discussion, and we set X = from this 
point onward. 

Given the solution for P{t^V), we can integrate over the bound states to find the fraction of 
systems in resonance as a function of time, i.e., 
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where k is the value of V for which the oscillator begins circulating. When k ^ Yt, the distribution 
function is primarily determined by diffusion, and the bound fraction takes the approximate form 

k 



Ph{t) = Erf 



(90) 



.2VZt, 

as in the stochastic pendulum problem. However, when Yt ^ k, at late times we can approximate 
the error function with the asymptotic formula. 
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which is vahd in the Umit x ^ 1 (Abramowitz & Stegun 1970). In this hmit, the fraction of bound 
states reduces to the form 



1 1 



f)"^7l7.-p(-5*). (.2) 



Yt-k Yt + k 

The fraction of bound states thus displays exponential decay in this hmit. 

In the limit where 62 <C e, which holds for the GJ876 system, the constants Y and Z reduce 
to the forms 



Y = -Pe. 7^ TT^ and Z = V, \ ] , (93) 



and hence the long term time evolution of the bound fraction is given by 



Pb ~ exp 



\{e) 



(94) 



From the definition of (5e2, we see that (5e2/(e)2 will be small when /Csi is small, i.e., when the 
system is not highly interactive. 

We can now summarize the implications of this solution: If the system is not sufHciently 
interactive, then the eccentricities do not vary appreciably; in mathematical terms, the factors 
5e2 and 5e will be small, the parameter Y will be small, and the system will exhibit a power-law 
decrease in the fraction of bound states as implied by equation ()90p . In the opposite limit of highly 
interactive systems, 5e2 and 5e are nontrivial, y 7^ 0, and the solution displays the exponential 
decay indicated by equation (j94p in the limit of long times. In this interactive case, since the 
eccentricities are varying and angular momentum is conserved, the semimajor axes of the planets 
will vary. As a result, the mean motions (17, Q2) will also vary and the location of the resonance 
will move around. This movement of the resonance thus corresponds to the "drift" in V obtained 
in the solution to the Fokker-Planck equation in the above analysis. When, in addition, the planets 
have large enough masses so that the star itself moves substantially (as in the case of GJ876), the 
location of the resonance moves even more. 

Next we want to compare the time scales for diffusion and drift. In the absence of the drifting 
effect, equation (j90p shows that the bound fraction evolves on a time scale idiffuse = A;^/(4Z). 
Furthermore, since Pb ~ in order for turbulence to have an effect, the evolution time must 

take place over many diffusion times; for example, Pb ~ 0.1 requires that t ~ lOOtdiffuse- When the 
drift term becomes significant, equation (f92l) shows that the bound fraction decays exponentially 
with a time scale texp = 4Z/y^. To simplify this discussion, let e and /x represent the eccentricities 
and mass ratios of both planets, and let 5 = 5ej/{e)j be the amplitude for eccentricity variations 
for both planets. The ratio of the two time scales is then given by 



iexp V 4^ / fJ- 
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Fig. 6. — Comparison between numerical integrations of the GJ876 system and the analytic results 
of this section. The solid curve shows the numerical results for the fraction of bound states, i.e., the 
number of planets in resonance divided by the number of planets in orbit. The level of turbulence 
is given by {Av)/v = 0.0002. The dashed curve shows the result expected from equation ([89]) . 
where the parameters are adjusted to provide a good fit. Note that the figure is presented as a 
log- linear plot, so that a straight line corresponds to exponential decay. The error bars show the 
root-A^ fluctuation amplitudes for comparison. 
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where we have absorbed numerical factors of order unity into the constant B. For relatively large 
eccentricities, say, e ~ 0.3 near the peak of the observed distribution for extrasolar planets, the 
factor (e^/fJ-) is of order unity (recall also that e = 0.26 for GJ876). If the diffusion constants 
are comparable, then equation ([95]) implies tdiflfuse/^exp ~ 5'^- Thus, in order for exponential decay 
to be realized, the parameter 5 must be relatively large, i.e., the changes in eccentricity must be 
comparable to the mean values (e). This condition is met when the two planets are relatively 
large and have comparable mass, so that Cs2 ~ ^ ~ ^/^ (^^^ equations [79] and [80]). 

Because exponential decay is much faster than power-law decay, this requirement on 6 is less severe 
then it might seem: Suppose, for example, turbulence acts for 100 diffusion times, so that the 
bound fraction Pb ~ 0.1 in the absence of interactions. The decay term is then roughly given by 
exp[— lOOtdiffuse/^exp] ~ exp[— 1005^]. Thus, the fraction of bound states can be affected for more 
moderately interacting systems, e.g., with values of 5 > 0.1. Keep in mind, however, that the ratio 
of time scales depends on the other variables (e, /x, !D„, Pet;) and will thus vary greatly from system 
to system. 

To test the predictions of this model, we compare the fraction of bound states (resonant 
systems) from equation (j92p with the fraction derived from an ensemble of numerical integrations 
of the GJ876 system. The result is shown in Figure 6, which shows the bound fraction - the ratio 
of the number of systems in resonance to the number of systems left in orbit - as a function of 
time. Both the analytic model and the numerical integrations show a nearly exponential decline 
in the number of bound states (Figure 6 is presented as a log-linear plot so that exponential 
decay corresponds to a straight line). The analytic model is given in terms of error functions with 
dimensionless arguments of the form x = {k i: Yt)/2^/Zt. As a result, any fit to the numerical 
result corresponds to a one parameter family of values of the parameters {k,Y,Z); the values will 
also depend on the choice of units. For the sake of definiteness, we have chosen k = 1. Since k is 
approximately the libration frequency of the resonance, the other two variables, here the diffusion 
constants Y and Z are then given in terms of libration frequencies. In these units (and for time t 
measured in years) we find that the values Y ~ 5.3 x 10~^ and Z 3.1 x 10~'^ provide a good fit 
(as shown in Figure 6). 

These results are sensible, at least in order of magnitude: From the original diffusion equation, 
the diffusion constant Z ~ ^o/t, where loq is the libration frequency and t is the diffusion time. 
The original stochastic differential equation also defines the diffusion time, i.e., the time required 
for random perturbations of size 77^ to change the libration speed by a total displacement of order 
wq. This time scale is given by t ~ 9Pu;q^I~'^[{Av) / v]'^ , where P is the period of the inner planet 
(in years). Putting these two results together, and inserting numerical values, we find that the 
diffusion constant is given by Z ~ 4(9^)^"*^ x 10~'^ ~ 5 x 10^'^ Y^^^- Note that this value must 
be adjusted to units in which A; = 1 to compare with the fitting values of our analytic theory to 
the numerical data; here, the libration period of GJ876 is about 9 years, so that k = 2-k/P ~ 0.7, 
which is close to unity, so that order of magnitude agreement obtains. 

One key result emerging from this analysis is that turbulence leads to a power-law decrease in 
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t (yr) 

Fig. 7. — Results from numerical integrations of a less interactive version of the GJ876 system 
where the inner planet has 10 times smaller mass than the observed system. The top panel shows 
the fraction of bound resonant states as a function of time during the early phase, when interactions 
have not yet had time to act and the bound fraction ~ t~'^l'^ . The dotted line shows a power-law 
fit to this result (the top panel is a log-log plot). The bottom panel shows the longer term trend (as 
a log-linear plot), when the bound fraction decreases exponentially with time; the dotted line shows 
a fit to this result. The dashed curve (bottom panel) shows the result expected from equation (j89p . 
where the parameters are adjusted to provide a good fit; the results of this section thus explain 
both the early power-law behavior and the later exponential behavior with a single model. 
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the fraction of resonant bound states in the regime of little interaction between planets, whereas 
highly interactive systems display exponential decay. As a test of these ideas, we have run the 
following numerical experiment: We begin with an ensemble of GJ876 system analogs, but reduce 
the mass of the inner planet by a factor of 10 in order to reduce the level of planet-planet interactions. 
All of the other system parameters are kept the same as before, and all of the systems are started 
in the 2:1 mean motion resonance. The level of turbulence, as set by the amplitude {l^v)/v of 
velocity perturbations, is also the same as before. 

The resulting time evolution for the fraction of bound states is shown in Figure 7 for an 
ensemble of A^ens = 10^ systems. The top panel shows the early time evolution as a log-log plot, so 
that power-law decay corresponds to a straight line. During this early phase of evolution, planetary 
interactions (which are much weaker than in the true GJ876 system) do not yet have time to act 
and cannot lead to drift of the resonance location. As a result, the straight dotted line in the top 
panel provides a good fit to the numerical results, with the expected power-law slope of —1/2 (see 
§2 and Paper I). At later times, the effects of planet-planet interactions gradually accumulate, and 
their effect on the dynamics becomes important. As shown above, at this stage, the time evolution 
of the fraction of bound states turns over to an exponentially decreasing form. This behavior is 
illustrated in the bottom panel of Figure 7, which is presented as a log-linear plot, so that straight 
lines correspond to exponential decay. Here, the dotted line shows a purely exponential decay, 
which is expected for the late time behavior of these systems. The dashed curve in the bottom 
panel shows the solution from equation ([89]) . where the parameters are adjusted to provide a good 
fit. Note that the model developed herein smoothly connects the non-interactive regime, with a 
power-law decrease in the fraction of resonant states, with the highly interactive regime, with an 
exponential decrease in the fraction of bound states. 

In order to derive the simplified Fokker-Planck equation (j87p and its solution for the time 
evolution of the fraction of bound states (equation [89]), we have made a number of approximations, 
and it is useful to summarize them here: First, we averaged over the libration angle (j) in the Fokker- 
Planck equation; this approximation is well known and well tested (e.g., MM04 and Paper I). At the 
next stage of the calculation, we adopted the ansatz of equations (I79p and ()80p . This approximation 
captures the basic behavior of the eccentricity variations for purposes of modeling the diffusion of the 
energy of the resonance; if we needed the full time dependence of the eccentricities, however, their 
full equations of motion should be retained. With this simplified description for the eccentricities, 
we then averaged over the apsidal angle 9 = V02 — vj, and thereby obtained a diffusion equation 
in two variables (e and V). Although the resulting two dimensional diffusion equation can be 
solved, it is rather cumbersome. In addition, the eccentricity values are not free to diffuse to 
arbitrary values, but rather are confined to the range < e < emax; where planets with e > Cmax 
tend to experience orbit crossing and ejection. We thus average the diffusion equation over the 
eccentricity to obtain the simplified form of equation ()86p . Although this procedure represents an 
uncontrolled approximation, and hence this last averaging is the least justified, the uncertainties 
can be encapsulated in the constants A and B. The solutions to the resulting diffusion equation can 
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then be found exactly, and they contain additional behavior, beyond that of the solutions found in 
§2. In particular, these solutions show that the fraction of bound states can decrease exponentially 
when the interactions between planets are significant, and this trend agrees with the results of the 
corresponding numerical simulations of the full 3-body problem. Finally, we note that equation 
(|87p is more robust than its derivation - this form of the Fokker-Planck equation represents a 
straightforward generalization of that obtained earlier for the simpler stochastic pendulum. 



5. CONCLUSION 

5.1. Summary of Results 

Building on earlier work (Paper I), this paper explores the effects of stochastic fluctuations on 
the maintenance of mean motion resonance in planetary systems. The principal result of this study 
is that turbulence can readily compromise mean motion resonance for the fluctuation amplitudes 
predicted by MHD simulations and for the solar system architectures observed in extrasolar plane- 
tary systems. Furthermore, we can understand the physical processes involved using a (primarily) 
analytic approach. A more specific outline of our results is given below: 

The most important quantity calculated in this paper is the fraction of systems that remain 
in resonance as a function of time. This fraction of bound states Ph{t) is defined to be the number 
of systems in resonance divided by the number of systems that remain intact (without ejecting 
a planet). This distinction is necessary because mean motion resonance protects multiple planet 
systems from instability, so that planetary ejection takes place with significant probability after a 
mean motion resonance is compromised. As a result, the fraction of systems remaining in resonance 
depends on whether the fraction is calculated relative to the original number of systems in the 
ensemble or the number of multiple planet systems remaining (Figure 1). 

The effects of stochastic fluctuations on mean motion resonance act in qualitatively different 
ways in systems that are highly interactive and those which are not (compare with Figure 2 with 
Figure 3). In systems with relatively little interaction between the planets, the fraction of bound 
states (systems in resonance) decreases with time as a power-law, specifically P\, ~ Systems 
that display this type of behavior are close to the idealized, circular restricted 3-body problem; the 
example considered in this paper contains a Jupiter in a nearly circular orbit (initially) in resonance 
with a Super-Earth on an interior orbit. The simple pendulum model of mean motion resonance 
(MD99) with turbulent forcing (Paper I) also derives from the circular restricted 3-body problem 
and leads to this same power-law time dependence (see equation [IB]). For these non- interactive 
systems, ensembles of the stochastic pendulum (§2.1 and 2.2), solutions to the Fokker-Planck 
equation (§2.3), and direct numerical integrations of the 3-body problem (§2.4) all predict the same 
time dependence for the fraction of systems remaining in mean motion resonance. In addition, 
the departures of the numerical results from the expectation values are well-characterized by the 
amplitude of root-A^ fluctuations (Figure 3). For this regime, we can summarize these results by 



- 36 - 



writing the expected fraction of surviving resonances in the form ^ C / N^rh ' , where Nq,.]-, is 
the total number of orbits for which turbulence is active, and where the dimensionless constant C 
depends on the amplitude of the fluctuations and the probability pret of returning to resonance after 
leaving (we expect C ~ 10 — 100; see Paper I). Thus, for A'^orb ~ 10^, we expect Pb ~ 0.01 — 0.1. 

In highly interactive systems, the fraction of bound states (systems in resonance) decreases 
more rapidly with time (see Figure 2) and eventually shows an exponential decay. In this paper, 
we have developed a generalized treatment for mean motion resonance that includes planetary 
interactions in the model equations (§4). This treatment initially retains six variables, and thus 
contains more physics than the simple one-variable pendulum model; in particular, we include 
terms that describe the interaction between planets due to their mutual excitation of eccentricities. 
In this model, the fraction of bound states Ph{t) shows two limiting regimes of behavior: For the 
case of minimal planet-planet interactions, or for sufficiently short time scales while the diffusion 
effects dominate, the bound fraction shows power-law behavior Pb ~ t"^^"^. For highly interactive 
systems, or for sufficiently late times, the evolution of the distribution function is dominated by 
drift terms due to interactions, and the model predicts an exponential decrease in the fraction of 
bound states (equation [S]). This predicted exponential behavior is in good agreement with that 
indicated by numerical experiments of highly interactive systems such as GJ876 (see Figures 2 and 
6). For these systems, the variations in the numerical results for the bound fraction Pb(^) are 
somewhat larger than expected from root-A^ fluctuations alone; this complication is most likely due 
to the planetary interactions, which allow for additional degrees of freedom. 

Systems that contain moderate levels of planet-planet interactions can fall in an intermediate 
regime, where the fraction of bound states initially decreases as a power-law with Pb ~ t^^/^ for 
a substantial time interval, but eventually switches over to an exponential decay (see Figure 7). 
Furthermore, the duration of the initial power-law phase depends on the level of interactions in 
the system. During the early time evolution, interactions have little effect, and the system acts 
essentially like the simple stochastic pendulum; at later times the effects of interactions accumu- 
late, and the fraction of bound states decays exponentially. These two regimes are connected at 
intermediate times, when Pb(i) behaves as a steeper power-law in time (see equation [92]). The 
model developed in §4 accounts for this early power-law behavior, the later exponential decay, and 
the smooth matching between the regimes. 

Finally, for completeness, we have also studied the ramifications of including a damping term 
in the model equations for mean motion resonance (§3). For example, this damping can be driven 
by torques from the circumstellar disk that is responsible for planetary migration. In this case, for 
sufficiently short time scales, the distribution function for the ensemble of states evolves like that of 
the stochastic pendulum; the energy expectation value increases linearly with time (equation |32j ) 
and the fraction of bound states decreases as a power-law Pb ~ t^^/^ (equation [39]). At later times 
when jt ^ 1, however, both the energy expectation value (equation [32]) and the fraction of bound 
states (equation [lO]) approach asymptotic values. The analytic solutions to the Fokker-Planck 
equation (§3) are in excellent agreement with results obtained from numerical integrations of the 
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stochastic pendulum equation with damping (see Figures 4 and 5). 

5.2. Discussion 

This paper extends the analysis of Paper I and bolsters its main conclusion, i.e., that turbulence 
implies mean motion resonances in extrasolar planetary systems should be relatively rare (roughly 
at the level of a few percent), unless turbulence has a limited duty cycle. The fact that turbulence 
is capable of removing systems from resonance is not surprising. As shown in equation ([TD, the 
potential energy associated with a bound resonant state is much smaller (typically, by a factor 
of ~ 10^) than the binding energy of the planet within the gravitational potential well of the 
star. Since planets can be ejected with relative efficiency during the early phases of solar system 
evolution (e.g., Rasio & Ford 1996), it makes sense that additional perturbations (here, turbulent 
fluctuations) can remove systems from resonance. 

Given that the results of this paper (and Paper I) suggest that mean motion resonances should 
be rare, it is useful to compare this prediction to existing observational data. However, we first note 
that a detailed comparison is premature, due to the small number of systems found to date, and 
due to selection effects. In the current sample, 30 systems are observed to have multiple planets. 
As outlined in the introduction, one system (GJ876) is found to be deep in a 2:1 resonance, whereas 
three other systems are either in resonance or close. If all three of these latter systems are actually 
in resonance, one "estimate" for the nominal rate of resonances would be 4/30 or 13%; if only 
GJ876 is truly in resonance, this rate becomes 1/30 or 3%. These estimates suggest that mean 
motion resonances are at least somewhat rare. Besides the small numbers involved, this percentage 
is uncertain for several reasons: One important issue is that the fraction of systems in mean motion 
resonance calculated here is the survival rate, i.e., the fraction of systems that start in resonance 
and are not removed via turbulence. A large fraction of the systems that leave resonance eject 
one of the planets, and hence would not be included in the number of observed multiple planet 
systems. Acting in the other direction, another consideration is that not all of the multiple planet 
systems in the current sample were ever (necessarily) in resonance. In addition, the observations 
are not complete, so that many of the single planet systems in the current sample might have 
companions (perhaps with nearly integer period ratios), which could either add to the number of 
resonant systems or add to the number that are not in resonance. Our understanding of these 
issues will benefit from future observations and hence better statistics. 

This work also shows that for systems that leave resonance, experience orbit crossing, and eject 
a planet, the surviving planet typically displays an eccentric orbit. Although a detailed statistical 
description of the resulting distribution is beyond the scope of this paper, we note that the full 
range of possible eccentricities is realized. This finding is consistent with previous simulations of 
multiple planet systems undergoing inward migration; in this setting, one of the planets is often 
scattered out of the system, and the remaining planets (for an ensemble of systems) are left with 
semimajor axes and eccentricities that fill the (a, e) plane in a manner consistent with the current 
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observational sample (Adams Sz Laughlin 2003, Moorhead &; Adams 2005). Preliminary numerical 

experiments of this process including turbulent fluctuations (Moorhead 2008) indicate that the 
orbital elements of the surviving planets continue to fill the (a, e) plane, but better statistics are 
needed (both from the simulations and for the observational sample). In any case, planetary 
scattering and ejection — perhaps enhanced by turbulence removing systems from resonance — 
provides one viable mechanism to increase the orbital eccentricity of extrasolar planets. 

Extrasolar planetary systems display a great deal of diversity, and the effects of turbulence 
will be different for varying solar system architectures. The results of this paper show that highly 
interactive systems (like GJ876) can reach a state of evolution where the fraction of resonant 
systems decays exponentially. Highly interactive systems are thus the least likely to be able to 
remain in resonance. This finding, however, makes the existence of GJ876 itself, which is observed 
to be deep in resonance, an even greater enigma. One way to account for the existence of the 2:1 
resonance in the GJ876 system is for its disk to have had low mass (to reduce the level of turbulent 
fluctuations) and/or a short survival time after planet formation (to reduce the duty cycle of the 
effect). However, this explanation for the survival of the resonance would make it more difficult for 
the disk to produce relatively large planets. Recall that both theoretical considerations (Laughlin et 
al. 2004) and observational findings (e.g., Johnson et al. 2007) suggest that M stars have difficulty 
forming planetary companions with Jovian masses. 

The results of §4 suggest that systems with the smallest levels of interaction between planets 
would have the best chance of surviving in a resonant state. However, this hypothesis is not entirely 
true: For a given orbital spacing, for example that corresponding to a 2:1 period ratio, the level 
of interactions decreases as the planetary masses decrease. On the other hand, planets with the 
smallest masses produce the smallest gaps in their circumstellar disks, and hence experience greater 
levels of turbulent forcing; in other words, the disk reduction factor F/j due to gaps is smaller for low 
mass planets (Paper I). Planetary systems with only lower mass planets, like the newly discovered 
HD 40307 system with masses of mp = 4 - 9 Me (Mayor et al. 2008), may experience greater 
levels of turbulent forcing than systems with larger planets, and hence would be more likely to be 
removed from resonance. Taken together, these theoretical results to date suggest that the systems 
with the best chances of maintaining mean motion resonance in the face of turbulence are those 
with planets of moderate mass and nearly circular orbits. The "Super-Earth" system of this paper 
— with a Jovian planet in an outer, nearly circular orbit and a smaller inner planet — provides 
one good example (see Figure 3). 

The numerical simulations indicate another difference between the highly interactive and less 
interactive systems. After leaving a resonant state, the Super-Earth systems take a long time for 
their planets to experience orbit crossing and eventual planetary ejection. The GJ876-type systems 
eject their planets much more readily. As a result, the less interactive systems (e.g., Super-Earth 
systems) can stay "close" to resonance for a long time, where the period ratios are close to 2:1, but 
the other orbital elements do not have the proper values to be in true resonance. 
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This paper has made progress on our understanding of how well the simplest pendulum equa- 
tion works as a model for mean motion resonances, including stochastic forcing terms, and we 
have generalized the model to include damping (§3) and planetary interactions in an approximate 
manner (§4). Nonetheless, a great deal of work remains to be done. One important issue is to 
understand how turbulence acts in the earliest stages of planet migration, i.e., when the planets are 
formed and first become locked into resonance. This issue requires much more extensive numerical 
work and is thus beyond the scope of this paper. Recent work (sec Moorhead 2008) is starting 
to study how turbulence affects the early phases of migration. The issue of turbulence during 
planetary formation is also being considered elsewhere (e.g., Masahiro et al. 2007). In addition, 
we have not included the back reaction of how planets (e.g., through gap clearing) can affect the 
generation of turbulent fluctuations through the magneto-rotational instability. All of these issues 
provide fruitful avenues for future research. 

In the coming years, as the statistics of observed multiple planet systems become sufficiently 
complete, the fraction of systems in resonance will provide important constraints on their formation 
and subsequent evolution. Since mean motion resonance can be compromised relatively easily, 
as shown herein, planetary systems observed in such bound states must have followed restricted 
historical paths (e.g., with little interaction and/or large dissipation). This type of analysis has 
been considered for the resonances in our Solar System (Peale 1976, Malhotra 1993) and for the 
GJ876 system (e.g., Lee & Peale 2002), and will soon provide interesting constraints on other 
extrasolar planetary systems. However, the greatest wealth of information on this topic will come 
when the observed sample of multiple planet systems is large enough to determine the fraction of 
actual mean motion resonances and the fraction of near-resonances. 
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